library(mediation)
library(foreign)
setwd("~/Replication Code/")

##---Loading in data ----
ccap0812 <- read.dta("./data/ccap0812_subset.dta")

vsg.dat <- read.dta("./data/vsg_subset.dta")

set.seed(1693)

#----------------------------------------------------------------------------------------------------#
# CCAP 2008-2012
#----------------------------------------------------------------------------------------------------#
dat1 <- na.omit(ccap0812[c("weight", "rr_12_sc", "rr_m_08_sc",
                            "obama_fav_m",
                            "pid7r_12_sc", "pid7r_m_08_sc")])

# Explaining 2008 politician evals
m.bo.ccap <- lm(obama_fav_m ~ pid7r_m_08_sc + rr_m_08_sc, 
                dat1, weights = weight)
summary(m.bo.ccap)

# PID on RR
m1.ccap <- lm(rr_12_sc ~ pid7r_m_08_sc + rr_m_08_sc, 
              dat1, weights = weight)
summary(m1.ccap)
m1.bo.ccap <- lm(rr_12_sc ~ obama_fav_m + pid7r_m_08_sc + rr_m_08_sc, 
                 dat1, weights = weight)
summary(m1.bo.ccap)

## ACME Mediation 
med.pid.ccap <- mediate(m.bo.ccap, m1.bo.ccap, 
                        treat = "pid7r_m_08_sc", 
                        mediator = "obama_fav_m",
                        boot = T, 
                        sims = 2000)
summary(med.pid.ccap)
sens.pid.ccap <- medsens(med.pid.ccap, rho.by = 0.01)
summary(sens.pid.ccap)

# RR on PID
m10.ccap <- lm(pid7r_12_sc ~ pid7r_m_08_sc + rr_m_08_sc, 
               dat1, weights = weight)
summary(m10.ccap)
m10.bo.ccap <- lm(pid7r_12_sc ~ obama_fav_m + pid7r_m_08_sc + rr_m_08_sc, 
                  dat1, weights = weight)
summary(m10.bo.ccap)


## ACME Mediation 
med.rr.ccap <- mediate(m.bo.ccap, m10.bo.ccap, 
                       treat = "rr_m_08_sc", 
                       mediator = "obama_fav_m",
                       boot = T, 
                       sims = 2000)
summary(med.rr.ccap)
sens.rr.ccap <- medsens(med.rr.ccap, rho.by = 0.01)
summary(sens.rr.ccap)

#----------------------------------------------------------------------------------------------------#
# VSG 2012-2016
#----------------------------------------------------------------------------------------------------#
dat2 <- na.omit(vsg.dat[c("rr_sc_12", "rr_sc_16", "pid7_sc_12", "pid7_sc_16", "weight",
                           "bo_fav_sc_12")]) 

# Explaining 2012 politician evals
m.bo.vsg <- lm(bo_fav_sc_12 ~ pid7_sc_12 + rr_sc_12, 
               data = dat2, weights = weight)
summary(m.bo.vsg)

# PID on RR
m1.vsg <- lm(rr_sc_16 ~ pid7_sc_12 + rr_sc_12, 
             data = dat2, weights = weight)
summary(m1.vsg)
m1.bo.vsg <- lm(rr_sc_16 ~ bo_fav_sc_12 + pid7_sc_12 + rr_sc_12, 
                data = dat2, weights = weight)
summary(m1.bo.vsg)

## ACME Mediation 
med.pid.vsg <- mediate(m.bo.vsg, m1.bo.vsg, 
                       treat = "pid7_sc_12", 
                       mediator = "bo_fav_sc_12",
                       boot = T, 
                       sims = 2000)
summary(med.pid.vsg)
sens.pid.vsg <- medsens(med.pid.vsg, rho.by = 0.01)
summary(sens.pid.vsg)

# PID on RR
m10.vsg <- lm(pid7_sc_16 ~ pid7_sc_12 + rr_sc_12, 
              data = dat2, weights = weight)
summary(m10.vsg)

m10.bo.vsg <- lm(pid7_sc_16 ~ bo_fav_sc_12 + pid7_sc_12 + rr_sc_12, 
                 data = dat2, weights = weight)
summary(m10.bo.vsg)

## ACME Mediation 
med.rr.vsg <- mediate(m.bo.vsg, m10.bo.vsg, 
                      treat = "rr_sc_12", 
                      mediator = "bo_fav_sc_12",
                      boot = T, 
                      sims = 2000)
summary(med.rr.vsg)
sens.rr.vsg <- medsens(med.rr.vsg, rho.by = 0.01)
summary(sens.rr.vsg)


#----------------------------------------------------------------------------------------------------#
# Plots
#----------------------------------------------------------------------------------------------------#
### Mediation effects
## PID on RR
med.pid <- array(NA, c(2,3))
med.pid.ci <- array(NA, c(2,6))
rownames(med.pid) <- c("CCAP", "VSG")
rownames(med.pid.ci) <- c("CCAP", "VSG")
colnames(med.pid) <- c("Total", "ACME", "ADE")
colnames(med.pid.ci) <- c("Total", "", "ACME", "", "ADE", "")
# Point Estimates
med.pid[1,3] <- med.pid.ccap$tau.coef
med.pid[2,3] <- med.pid.vsg$tau.coef
med.pid[1,2] <- med.pid.ccap$d.avg
med.pid[2,2] <- med.pid.vsg$d.avg
med.pid[1,1] <- med.pid.ccap$z.avg
med.pid[2,1] <- med.pid.vsg$z.avg
## CIs
med.pid.ci[1,5:6] <- med.pid.ccap$tau.ci
med.pid.ci[2,5:6] <- med.pid.vsg$tau.ci
med.pid.ci[1,3:4] <- med.pid.ccap$d.avg.ci
med.pid.ci[2,3:4] <- med.pid.vsg$d.avg.ci
med.pid.ci[1,1:2] <- med.pid.ccap$z.avg.ci
med.pid.ci[2,1:2] <- med.pid.vsg$z.avg.ci

## RR on PID
med.rr <- array(NA, c(2,3))
med.rr.ci <- array(NA, c(2,6))
rownames(med.rr) <- c("CCAP", "VSG")
rownames(med.rr.ci) <- c("CCAP", "VSG")
colnames(med.rr) <- c("Total", "ACME", "ADE")
colnames(med.rr.ci) <- c("Total", "", "ACME", "", "ADE", "")

# Point Estimates
med.rr[1,3] <- med.rr.ccap$tau.coef
med.rr[2,3] <- med.rr.vsg$tau.coef
med.rr[1,2] <- med.rr.ccap$d.avg
med.rr[2,2] <- med.rr.vsg$d.avg
med.rr[1,1] <- med.rr.ccap$z.avg
med.rr[2,1] <- med.rr.vsg$z.avg
## CIs
med.rr.ci[1,5:6] <- med.rr.ccap$tau.ci
med.rr.ci[2,5:6] <- med.rr.vsg$tau.ci
med.rr.ci[1,3:4] <- med.rr.ccap$d.avg.ci
med.rr.ci[2,3:4] <- med.rr.vsg$d.avg.ci
med.rr.ci[1,1:2] <- med.rr.ccap$z.avg.ci
med.rr.ci[2,1:2] <- med.rr.vsg$z.avg.ci

# Figure 1
# pdf("./figures/med_plots.pdf", width = 10, height = 5)
par(mfcol = c(1,2))
plot(med.pid[1,], 1:3,
     xlim = c(0,
              max(med.pid.ci) + .02),
     ylim = c(1,3.2),
     pch = 15,
     main = "Partisanship on Racial Attitudes",
     xlab = "Estimate",
     ylab = "",
     axes = F)
points(med.pid[2,], seq(1.1,3.1,1), pch = 16, col = "grey")
segments(y0 = 1:3,
         x0 = med.pid.ci[1,c(1,3,5)],
         x1 = med.pid.ci[1,c(2,4,6)])
segments(y0 = seq(1.1,3.1,1),
         x0 = med.pid.ci[2,c(1,3,5)],
         x1 = med.pid.ci[2,c(2,4,6)],
         col = "grey")
axis(1, at = seq(0, (max(med.pid.ci) + .02), .05))
axis(2, at = seq(1,3), labels = c("ADE", "ACME", "Total\nEffect"), las = 1)
legend("topleft", c("2008-12\nCCAP", "2012-2016\nVOTER Survey"),
       pch = c(15,16), col = c("black", "grey"), bty = "n",
       y.intersp = 2)

plot(med.rr[1,], 1:3,
     xlim = c(0,
              max(med.rr.ci) + .02),
     ylim = c(1,3.2),
     pch = 15,
     main = "Racial Attitudes on Partisanship",
     xlab = "Estimate",
     ylab = "",
     axes = F)
points(med.rr[2,], seq(1.1,3.1,1), pch = 16, col = "grey")
segments(y0 = 1:3,
         x0 = med.rr.ci[1,c(1,3,5)],
         x1 = med.rr.ci[1,c(2,4,6)])
segments(y0 = seq(1.1,3.1,1),
         x0 = med.rr.ci[2,c(1,3,5)],
         x1 = med.rr.ci[2,c(2,4,6)],
         col = "grey")
axis(1, at = seq(0, max(med.rr.ci) + .02, .05))
axis(2, at = seq(1,3), labels = c("ADE", "ACME", "Total\nEffect"), las = 1)
legend("topleft", c("2008-12\nCCAP", "2012-2016\nVOTER Survey"),
       pch = c(15,16), col = c("black", "grey"), bty = "n",
       y.intersp = 2)
# dev.off()

### Sensitivity Analysis (Figure 2)
# pdf("./figures/med_sens.pdf")
par(mar = c(5,6,3,1), mfrow = c(2,2))
plot(sens.pid.ccap, sens.par = "R2",
     r.type = 2, sign.prod = -1,
     xlim = c(0, 0.6), ylim = c(0, 0.4),
     main = "Partisanship's Mediated Effect\non Racial Attitudes\n2008-12 CCAP",
     xlab = "Proportion of Variance in Obama Evaluations\n Explained by Confounder",
     ylab = "Proportion of Variance in Racial Attitudes\n Explained by Confounder",
     axes = F,
     cex.main = .95) # sign.prod = 1 if confounding is same direction for mediator and outcome; -1 if different
plot(sens.rr.ccap, sens.par = "R2", 
     r.type = 2, sign.prod = -1,
     xlim = c(0, 0.6), ylim = c(0, 0.3),
     main = "Racial Attitudes' Mediated Effect\non Partisanship\n2008-12 CCAP",
     xlab = "Proportion of Variance in Obama Evaluations\n Explained by Confounder",
     ylab = "Proportion of Variance in Partisanship\n Explained by Confounder",
     axes = F,
     cex.main = .95)
plot(sens.pid.vsg, sens.par = "R2",
     r.type = 2, sign.prod = -1,
     xlim = c(0, 0.5), ylim = c(0, 0.4),
     main = "Partisanship's Mediated Effect\non Racial Attitudes\n2012-16 VOTER Survey",
     xlab = "Proportion of Variance in Obama Evaluations\n Explained by Confounder",
     ylab = "Proportion of Variance in Racial Attitudes\n Explained by Confounder",
     axes = F,
     cex.main = .95) 
plot(sens.rr.vsg, sens.par = "R2", 
     r.type = 2, sign.prod = -1,
     xlim = c(0, 0.5), ylim = c(0, 0.3),
     main = "Racial Attitudes' Mediated Effect\non Partisanship\n2012-16 VOTER Survey",
     xlab = "Proportion of Variance in Obama Evaluations\n Explained by Confounder",
     ylab = "Proportion of Variance in Partisanship\n Explained by Confounder",
     axes = F,
     cex.main = .95)
# dev.off()